function   imp = irf_to_dN(    current_model , years_simulation, w, lambda,  options )

if( nargin<=5 )
    options.random = false;
end

disp(  [  '----------------- Generating the impulse response function for a ', current_model.par.model_name, ' model --------------------'   ]  )

par = current_model.par;
% Find the correct scaling for the credit spread. 
w_dim = kron(ones(par.N_lambda,1), par.w_grid );
lambda_dim = kron(par.lambda_grid, ones(par.N_w,1));
fit_a_function = @(values_matrix)  scatteredInterpolant(w_dim,    lambda_dim ,  reshape( values_matrix, par.N_w*par.N_lambda, 1 ),   'linear')  ;
m = simulation_results_output(current_model);
spread_fitted = fit_a_function(current_model.credit_spread);
avg_spread_before_normalization = mean( spread_fitted(m.w,m.lambda) ) ;
spread_fitted = fit_a_function(current_model.credit_spread/avg_spread_before_normalization);   % Re-adjust with the average

dt=par.dt;
N_sim = round(years_simulation/dt);

dN_vec = zeros( N_sim, 1); 
dN_vec(1) = 1;

% Then start the simulation. 
if( ~options.random )  % If we do not want shocks after the simulation starts    
    dB_vec = zeros( N_sim, 1 ); 

    
    sim = simulate( w, lambda, dB_vec, 0*dN_vec,  years_simulation, current_model,  par  );
    sim_pre_recap = sim.results;
    
    sim = simulate( w , lambda, dB_vec, dN_vec,  years_simulation, current_model,  par  );
    sim_post_recap = sim.results;
    
    imp.w_pre = sim_pre_recap.w;   imp.w_post=sim_post_recap.w;
    % Impulse responses
    w_impulse = sim_post_recap.w ./ sim_pre_recap.w - 1;
    cs_impulse =  spread_fitted( sim_post_recap.w,   sim_post_recap.lambda ) - spread_fitted( sim_pre_recap.w,   sim_pre_recap.lambda );
    output_impulse =   sim_post_recap.GDP ./ sim_pre_recap.GDP -1  ;
    bank_credit_impulse = sim_post_recap.bank_credit ./ sim_pre_recap.bank_credit - 1;
else    
    if( isfield(options, 'N_rounds' ) )
        N_rounds = options.N_rounds;    
    else
        N_rounds = 400;    
    end
    
    w_impulse_matrix = zeros(N_sim, N_rounds);   cs_impulse_matrix = zeros(N_sim, N_rounds); 
    output_impulse_matrix = zeros(N_sim, N_rounds);   bank_credit_impulse_matrix = zeros(N_sim, N_rounds);
    output_impulse0_matrix= zeros(N_sim, N_rounds);
    
    parfor( iter=1:N_rounds  )
        if( mod( iter, 50 )==0 )
            disp( ['Progress: ', num2str(iter/N_rounds*100), '%']  )
        end        
        results = shocks_generation( years_simulation,  par,  iter*10 );
        dB_vec = results.dB_vec;

        sim = simulate( w, lambda, dB_vec, 0*dN_vec,  years_simulation, current_model,  par  );
        sim_pre_recap = sim.results;

        sim = simulate( w , lambda, dB_vec, dN_vec,  years_simulation, current_model,  par  );
        sim_post_recap = sim.results;
        
        % Impulse responses
        w_imp = sim_post_recap.w ./ sim_pre_recap.w - 1;
        cs_imp =  spread_fitted( sim_post_recap.w,   sim_post_recap.lambda ) -  spread_fitted( sim_pre_recap.w,   sim_pre_recap.lambda ) ;
        if( max(max(current_model.credit_spread/avg_spread_before_normalization))>2 )
            cs_imp = cs_imp/2;
        end
        if( max(cs_imp)>0 )
            disp( ['credit spread response abnormal at iter=', num2str(iter)]  )
        end
        
        output_imp =   sim_post_recap.GDP ./ sim_pre_recap.GDP -1  ;
        output_imp0 = sim_post_recap.GDP ./ sim_pre_recap.GDP(1) -1  ;

        bank_credit_imp = sim_post_recap.bank_credit ./ sim_pre_recap.bank_credit - 1;
        
        w_impulse_matrix(:,iter) = w_imp;
        cs_impulse_matrix(:,iter) = cs_imp;
        output_impulse_matrix(:,iter) = output_imp;
        output_impulse0_matrix(:,iter) = output_imp0;
        bank_credit_impulse_matrix(:,iter) = bank_credit_imp;
    end
    w_impulse = mean(w_impulse_matrix,2) ;
    cs_impulse = mean(cs_impulse_matrix,2);
    output_impulse = mean(output_impulse_matrix,2);
    bank_credit_impulse = mean(bank_credit_impulse_matrix,2);
    
    imp.w_impulse_matrix = w_impulse_matrix;
    imp.cs_impulse_matrix = cs_impulse_matrix;
    imp.output_impulse_matrix = output_impulse_matrix;  
    imp.output_impulse0_matrix = output_impulse0_matrix;
    imp.bank_credit_impulse_matrix = bank_credit_impulse_matrix;
end

if(   isfield( options, 'smooth' ) && options.smooth  )
    x = [1:N_sim]*dt;
    process_fun = @(y) slmeval(  x,  slmengine( x, smooth(y) , 'plot', 'off' )  );
else
    process_fun = @(y)  y;
end

imp.w_impulse = process_fun(w_impulse);
imp.cs_impulse = process_fun(cs_impulse);
imp.output_impulse = process_fun(output_impulse);
imp.bank_credit_impulse = process_fun(bank_credit_impulse);
imp.lambda_init = lambda;



